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ABSTRACT 

CoUisionless simulations of the CDM cosmology predict a plethora of dark matter sub- 
structures in the halos of Milky Way sized galaxies, yet the number of known luminous 
satellites galaxies is very much smaller, a discrepancy that has become known as the 
'missing satellite problem'. The most massive substructures have been shown to be 
plausibly the hosts of the brightest satellites, but it remains unclear which processes 
prevent star formation in the many other, purely dark substructures. We use high- 
resolution hydrodynamic simulations of the formation of Milky Way sized galaxies in 
order to test how well such self-consistent models of structure formation match the 
observed properties of the Galaxy's satellite population. For the first time, we include 
in such calculations feedback from cosmic rays injected into the star forming gas by 
supernovae as well as the energy input from supermassive black holes growing at the 
Milky Way's centre and its progenitor systems. We find that non-thermal particle pop- 
ulations quite strongly suppress the star formation efficiency of the smallest galaxies. 
In fact, our cosmic ray model is able to reproduce the observed faint-end of the satellite 
luminosity function, while models that include only the effects of cosmic reionization, 
or galactic winds, do significantly worse. Our simulated satellite population approxi- 
mately matches available kinematic data on the satellites and their observed spatial 
distribution. We conclude that a proper resolution of the missing satellite problem 
likely requires the inclusion of non-standard physics for regulating star formation in 
the smallest halos, and that cosmic reionization alone may not be sufficient. 
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1 INTRODUCTION 

The leading ACDM cosmology predicts tha t galaxies form 
hierarchically in a 'bottom up' fashion (e.g. IWhite fc ReesI 
ll978l : lLongaiilll999l ). where small perturbations in the dark 
matter density distribution collapse earlier than larger per- 
turbations, and low-mass halos grow by smooth accretion 
or mergers with other halos, successively building up ever 
bigger structures. But structures falling into bigger systems 
during this process are not always disrupted completely. As 
N-body simulations show, the inner cores of infalling objects 
often survive the various disruptive effects acting on them, 
like tidal truncation, tidal shocking or ram-pressure strip- 
ping. It is believed that the observed dwarf galaxies orbiting 
around the Milky Way (MW) are examples of such surviving 
remnants. 

Based on the first generation of very high resolu- 
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tion coUision l ess C DM simulations, iKlvpin et al. I l|l999l ) and 
iMoore et al.1 (| 19991 ) pointed out a very striking apparent dis- 
crepancy between theoretical predictions for such satellite 
systems and actual observations. Given the very large num- 
ber of predicted dark matter substructures, there appears to 
be a dearth of luminous satellites in the Milky Way. In fact, 
the cumulative number of observed satellite galaxies and of 
predicted substructures above a given circular velocity value 
differed by a factor of ~ 10. This has become known as the 
'missing satellite problem'. 

The ini t ial a nalysis of iMoore et al.l (|l999l ) and 
iKlvpin et al.l (|l999l ) may have overstated the magnitude 
of the discrepancy, both because of uncertainties in as- 
signing co rrect circular velo city values to the observed 
satellites (jStoehr et al.l |2002| ) and because a number of 
additional faint satellites have been discovered mean- 
while in the M W (see for example 'irwin et al.' '2007; 



Liu et all ,2 008; Martin et al., 2008; Simon & Geha 2007|; 



2000i : Ivan den Berghl [2OO0I : iBelokurov et al.l 12003 . 
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I2OO9I : IWatkins' et al. l l2009l : lBelokurov et al.|[201ol '). However, 
there is a consensus that the many low-mass satelhtes pre- 
dicted by the N-body simulations need to be strongly sup- 
pressed in luminosity, otherwise a significant discrepancy 
with the observed satellite luminosity functions results, that, 
if confirmed, may in principle even be used to rule out cold 
dark matter. 

With increasing numerical resolution, the missing satel- 
lite problem has become more acute. Modern cosmolog- 
ical dark matter simulations of Milk y Wa y sized ha- 



plain the observed satellite abundance (|Maccid fc FontanotI 



los llDiemand et al.|[20o3 : ISpringel et al.N 2008': Sta del etal] 
I2OO9F resolve up to ^ 300, 000 dark matter substructures, 
while the number of observed satellite galaxies still com- 
prises just a few dozens. We note that the modern dark 
matter only simulations are even able to resolve substruc- 
tures inside substructures, and interestingly, there is also 
some observational evide nce for a satellite poss ibly orbiting 
aro und another satellite (|Belokurov et al.]|2008l ). 

iKlvpin et al.l (|l999l l did not only raise the missing satel- 
lite problem, they were also among the first to suggest 
a potential solution to this issue. In particular, they pro- 
posed that star formation inside low mass halos could be 
suppressed because of photo evaporation of gas due to a 
strong intergalactic ionizing UV background. This would 
keep most of the orbiting satellites dark and render them vi- 
su ally unobservab le. Indeed, the simple filtering mass model 
of iGnedinI (|2000l l for the impact of a UV background on 
the cooling efficiency of small halos predicts a quite sizable 
effect, with a nearly complete suppression of cooling in all 
halos with circular velocity below 50 km s~^ . However, recent 
full hydrodynamical simulations have not confirmed this 
l|Hoeft et al.ll2006l : [Okamoto et al.|[2008l ). They find a con- 
siderably weaker effect, where only halos with circular veloc- 
ities less than ^ 25kms~^ are affected. This also casts some 
doubt about the faint-end results of nu merous semi-analytic 
models for the satell ite population (e.g. iBenson et al.|[2002l : 
iKravtsov et al.ir2004f ) , which typically employed the filtering 
mass formalism and hence assumed an overly strong effect 
of the UV background. We will reexamine this question in 
this work based on our cosmological hydrodynamic simula- 
tions of Milky Way formation, which include a treatment of 
cosmic reionization. 

Another possible solut i on for the satellite problem was 
proposed by IStrigari et al.l (|2007l l who suggested that not 
the satellite mass at the present epoch determines whether a 
satellite would be luminous or not, but rather the maximum 
mass it had before accretion onto the Milky Way's halo. This 
is based on the idea that tidal stripping and ram pressure 
unbinds the gas from an infalling satellite and thus stalls 
any further star formation. With this assumption, the stellar 
mass at the time of accretion is essentially retained until the 
present epoch, and it becomes a question of allowing high- 
redshift star formation only in satellites above a sufficiently 
high mass threshold. 

A more radical conjecture is that the properties of the 
dark matter particles may have to be changed. Instead of 
having negligible velocity dispersion at the time of decou- 
pling, we may instead be dealing with (slightly) warm dark 
matter (WDM). This can supp ress the abundan ce of low 
mass structure considerably fe.g. lColm et al.ll2000l ). but pro- 
vided the particle mass is not lower than ~ 1 keV a suffi- 
ciently large number of substructures still survives to ex- 



In the most recent works on the subject, a num- 
ber of interesting and encouraging results have been ob- 
tained. Observationally, it has been discovered that the 
satellites all have approximately the same central mass den- 
sity (w ithin 300 to 600 p c ), independent of t heir lumi- 
nosity (|Gilmore et al.ll2007l : IStrigari et ai]|2008l ). Explain- 
ing this central density threshold has become an impor- 
tant additional challenge for theoretical models. Also, a 
significant number of new faint satellit es have been dis- 
covered with t he help of the SPSS Jlrwin et all l2007l : 



Liu et al. 2008 1: iMartin et all [200^ : ISimon fc Geha 
Grebeli ,200Q: Ivan den B ergh' ^200 0^ iBelokurov et al 



20091 : 



Watkins et al.l 2009; Bclokurov et alj 



2Q10|), 



2007; 



200a 



improv- 
ing our knowledge about the full satellite population sig- 
nificantly, but at the same time also raising the question 
whether we may perhaps still be missing large numbers of 
satellites at ultra low surface brightnesses. 

On the theoretical side, refined treatments of the effects 
of reionization, often coupled to the results of high-resolution 
coUisionless sim ulations have been u sed to model the satel- 
lite population. iMaccid et al.l (|2010l ) employed a number of 
different semi-analytic models and low-resolution hydrody- 
namic simulations to study the satellite luminosity function. 
Despite just invoking photoheating as primary feedback pro- 
cess, they achieved reasonable agreement for some of their 
models, leading the m to argu e that the satellite problem may 
be solved. Similarlv. lLi et al.l (|2010[ ) invoked a strong impact 
of reionization in a semi-analytic mo del similar to those ap- 
plied to the Millennium Simulation l|Croton et al.|[2006l ) to 
reproduce t he luminosity funct ion of galaxies around the 
Milky Way. iBusha et ail (|2010l ) used simple prescriptions 
for the impact of inhomogeneous reionization on the satellite 
population, pointing out that subtle changes in the assump- 
tions about how reionization affects star formation in small 
galaxies can lead to large changes in the predicted number 

of satellites. 

Recently, IStrigari et ai] (|2010l ') examined the kinematics 
of five well-measured Milky Way satellite galaxies and com- 
pared them to dark matter satellites of the high-resolutio n 
simulations of the Aquarius Project fepr ingel et al.|[2008h . 
They showed that these systems are fully consistent with 
ACDM expectations and may be hosted in cored dark mat- 
ter structures with maximum circular velocit i es in t he range 
10 to 30kms"\ Interestinglv. iBuUock et al.l (|2009l ') pointed 
out that the number of real satellite systems may in fact be 
much larger than commonly believed, with the majority of 
them being so far undetected because of their low surface 
brightness. In this scenario, the 'common mass scale' in- 
ferred for the observed satellites may in fact just arise from 
a selection bias. 

The first high-resolution hydrodynamic simulation able 
to directly r esolve the satellite popula tion has recently been 
presented bv lOkamoto fc FrenkI (|2009l '). They argue that the 
common mass scale identified in the observations arises from 
early reionization at redshift around z ~ 12, and that satel- 
lites that have not yet grown to a maximum circular velocity 
of ~ 12kms^^ by the time of reionization, will not be able 
to make any stars later on. Even if they grow above this 
threshold. . Okamoto fc Frenk ((200S ) predict them to remain 
dark. 
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Despite all of this progress, it is evident that there re- 
main many open questions concerning the population of 
faint and ultra-faint satellite galaxies orbiting around Milky 
Way like galaxies. Especially the influence of different bary- 
onic feedback processes on the luminosity function of the 
simulated satellites has not been investigated in sufBcient de- 
tail. It is therefore far from clear whether photoheating from 
a UV background and ordinary supernovae feedback are in- 
deed the correct physical solutions to the missing satellite 
problem. In fact, the in part contradictory results that have 
been obtained with analytic recipes to describe the impact of 
reionization suggest that more accurate methodologies are 
required to reliably settle the issue. 

We have therefore embarked on a research program 
where we use high-resolution hydrodynamic simulations of 
the formation of Milky Way-sized halos to shed more light 
on these questions, in particular by investigating a variety 
of feedback processes known to be important in galaxy for- 
mation. Besides the impact of reionization, these include 
galactic winds and outflows, energy input by growing super- 
massive black holes, or the non-thermal support of gas by 
cosmic rays or magnetic fields. Ultimately we aim to reach 
similar numerical resolution as has been obtained for recent 
collisionless simulations, even though this goal may still be 
several years away. 

In this work, we present some of our first results. We 
use several well resolved hydrodynamical simulations of the 
formation of a Milky Way sized galaxy to investigate the 
properties of the predicted population of satellite galaxies, 
for different choices of the included physics. Besides a default 
reference model that includes only a treatment of radiative 
cooling, star formation, and cosmic reionization, we consider 
also models that add galactic winds, supermassive black hole 
growth, or cosmic ray injection by supernovae shock waves. 
By comparing the simulation results with a comprehensive 
catalogue of the known Milky Way satellites, we seek to de- 
termine which of these processes is most important in shap- 
ing the satellite population. 

This paper is organized as follows. In Section [2] we de- 
scribe the methodological details of our simulations, while 
the observational knowledge about the satellites is briefly 
summarized in Section [S] Sections |4l [5] and [6] present the re- 
sults for our simulated populations of satellite galaxies, both 
with respect to individual satellite histories as well as with 
respect to their population as a whole. Our conclusions are 
summarized in Section [T] 



2 METHODOLOGY 

Our simulations are based on initial conditions originally 
constructed for the Aquarius Project (jSpringel et al.l 20081 ) 
of the Virgo Consortium. This project carried out highly 
resolved dark matter only simulations of 6 different Milky 
Way-sized halos, at a variet y of different numeric al resolu- 
tions. In the nomenclature of lSpringel et all (|2008l ). the halo 
and resolution level investigated here is called 'Aq-C-4', and 
has about 5.4 million dark matter particles in the flnal virial 
radius. The same object has also been stud i ed in the hydro- 
dynamic simulations of IScannapieco et al.l (|2009l ). albeit at 
the considerably lower resolution (by a factor 8 in particle 
number) corresponding to 'Aq-C-5'. There it was found that 



this 'C'-halo produced the lowest bulge-to-disk ratio among 
the 6 candidate halos selected in the Aquarius Project, mak- 
ing it a particularly good candidate for the formation of a 
large disk galaxy. We note, however, that we do not expect 
our choice of target halo to influence our principal conclu- 
sions for the satellite population. 

In this paper, we model the gas component with 
smoothed particle hydrodynamics (SPH) and introduce the 
gas particles into the initial conditions by splitting each orig- 
inal particle into a dark matter and gas particle pair, dis- 
placed slightly with respect to each other (at flxed center- 
of-mass) to arrive at a regular distribution of the mean par- 
ticle separations, and with a mass ratio corresponding to a 
baryon fraction of 16 per cent. The cosmological parameters 
are i},n ~ 0.25, Qa ~ 0.75, as = 0.9 and h = 0.73, the same 
ones used as in the original Aquarius simulations, which are 
consistent with the WMAPl cosmological constraints. A pe- 
riodic box of size 100/i"^Mpc on a side is simulated, with 
varying spatial resolution that 'zooms in' on the formation 
of a single galaxy. In the high-resolution region, we reach a 
mass simulation of ?a 2 x 10^ H'^Mq and 2 x lO"* H'^Mq 
for dark matter and gas particles, respectively. A constant 
comoving gravitational softening length of e = 0.25/i^^kpc 
was used for all high-resolution particles. 

We employed the parallel TreeSPH code gadget-3 for 
our runs, which is an imp roved and extended version of 
GADGET-2 l|Springell |2005| ). gadget calculates the long- 
range gravitational fleld in Fourier space, and the short 
range forces in real space with a hierarchical multipole 
expansion, based on a tree. This approach guarantees a 
homogeneously high spatial resolution in the gravitational 
force calculation and can be efficiently combined with an 
individual timestep integration scheme. For the hydrody- 
nam ics, gadget uses the 'entropy formulation' of SPH 
( Springel fc HernauistI |2002| ) , which is derived from a vari- 
ational principle and simultaneously conserves energy and 
entropy where appropriate. 

In the hydrodynamic part of gadget, different physi- 
cal processes besides ordinary gas dynamics are calculated. 
Most importantly, these are radiative cooling, star forma- 
tion and its regulation by supernovae feedback processes 
ISpringel fc HernauistI boOSH . The code can optio nally also 
model black hole growth and (ISpringel et al.ll2005l ) and cos- 



mic ray physics (jjubelgas et al. 2008l l . We shall now briefly 



describe the physics modules we used. 



2.1 Star formation model 

Radiative cooling is followed for a primordial mixture of 
helium and hydrogen under the assumption o f collisional 
ionization equilibrium, using a formulation as in iKatz et al.l 
(|l996l ). A spatially uniform, ionizing UV background is in- 
troduced with an amplit ude and time evolution described 
by an updated version of iHaardt fc Madaul (|l996l ). leading 
to reionization of the universe by redshift « ~ 6. 

To model star formation, we use the hybrid multiphase 
m odel for star formation a nd su pernova feedback introduced 
bv lSpringel fc HernauistI l|2003f ). in which every sufficiently 
dense gas particle is treated as a representative region for 
the multiphase structure of the interstellar medium (ISM). 
These hybrid particles are pictured to be comprised of cold 
dense clouds in pressure equilibrium with a hot ambient gas 
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phase, where only the clouds contribute material available 
for star formation. Mass and energy exchange processes be- 
tween these two phases a re computed by simple differ ential 
equations, as described in lSpringel fc IfernquistI l|2003l ). giv- 
ing rise to an effective equation of state that regulates the 
dense gas of the ISM. CoUisionless star particles are spawned 
stochastically from this star-forming phase according to a lo- 
cal estimate of the star formation rate. The gas consumption 
timescale of the m odel is calibrated such that it reproduces 
the Kennicutt law (|Kennicuttlll989l ) between star formation 
rate and gas surface density observed in low-redshift disk 
galaxies. 



2.2 Wind model 

Even though the above ISM model reproduces the star for- 
mation rates observed in disk galaxies, it overproduces the 
total amount of stars formed when applied in cosmological 
simulations. This is likely related to its inability to reproduce 
the observed galajcy-scale outflows seen around many star- 
bursting galaxies. We have therefore also investigated the 
phenomenological mo del for galactic winds introduced by 
ISpringel fc Hernguistl l|2005 ). We adopted a constant wind 
velocity of iiwind = 484kms~^ and mass loading factor of 
7] = 2, i.e. the mass flux of the wind is twice the star forma- 
tion rate. Individual gas particles were stochastically added 
to the galactic wind by changing their velocity to the pre- 
scribed wind velocity. We adopted an anisotropic distribu- 
tion for the wind direction, launching it preferentially per- 
pendicular to the disk. 

The wind model causes an outflow of gas from the dense 
gaseous disc, transporting energy, matter and heavy ele- 
ments out of the disc in proportion to the star formation 
rate. Not only the central galaxy is affected by the winds, 
all the star forming satellite galaxies produce winds, which 
will more quickly deplete their gas content. The effect of 
galactic winds is actually expected to be more effective in 
these low mass systems as their potenti al well is much sha l- 
lower than the one of the main galaxy (|Dekel fc Silklll986h . 
which should increase the mass loss due to outflows. In- 
cluding winds in the simulation thus appears in principle 
promising as a mechanism to reduce the number of lumi- 
nous satellites, but the effect may sensitively depend on 
how the wind properties a re scaled with galaxy size (e.g. 
lOppenheimer fc Dave|[200^ ). 



2.3 Black hole model 

Supermassive black holes are thought to reside at the cen- 
tres of most if not all spheroidal galaxies. The tight relation 
between their masses and the velocity dispersion of their 
hosting galaxies suggest a tight evolutionary link, which is 
probably established by a self-limited growth mechanism 
in which the energy output of a growing black hole even- 
tually terminates its further growth and the surrounding 
star formation, for example by expelling gas from the cen- 
tral region of the galaxy. Hydrodynamical simulation mod- 
els h ave been successfully used to model this pr ocess in de- 
tail (|Di Matteo et al.ll2005l : ISpringel et al.ll2005l 'l. and led to 
quite suc cessful unified model s of the formation of spheroidal 
galaxies l|Hopkins et al.ll2006l ). 



While it is unclear whether the influence of the Milky 
Way's black hole has affected the formation of other compo- 
nents of the Galaxy besides the central bulge, it appears pos- 
sible that the heating effects from different quasar episodes 
during the growth history of the MW's supermassive black 
hole have had an impact on the satellite population as well, 
for example by heating the environment of progenitor halos 
through strong quasar driven outflows. Indeed, in theoret- 
ical models outflow feedback has been found to be violent 
eno ugh to be able to strongly affect even neighbou ring galax- 
ies (IScannapieco et alllioOll : iThacker et al.il2002h . 

To study such effects, w e have adopted the techniques 
introduced by ISpringel et all {200^ for tracking black hole 
growth and its associated energy feedback in cosmological 
simulations. In brief, we periodically call a FoF group flnd- 
ing algorithm that identifles newly formed halos that do not 
contain a black hole yet. If such a halo is sufficiently mas- 
sive, its densest gas particle is converted to a black hole 
seed of mass Mbh = 10* h~^MQ. The black hole particles 
are treated as sinks particles that accrete gas from their 
surroundings at a rate estimate with a simple Bondi-Hoyle 
prescription, limited to the Eddington rate. The black hole 
accretion is assumed to have a fixed radiative efficiency of 
0.1 A/c^, and 5% of the produced radiation are assumed to 
couple thermally to the gas surrounding the black hole. This 
energy feedback can eventually drive a quasar-driven outflow 
and regulates the mass growth of the black holes. We also 
allow for two black holes to merge with each other once they 
get sufficiently close to each other. 



2.4 Cosmic ray model 

In the interstellar medium of our own Galaxy, it is believed 
that thermal pressure, cosmic rays and magnetic fields are 
roughly in equipartition. Even though it is hence known that 
non-thermal particle populations play an important role in 
regulating the gas dynamics of the ISM, this component has 
usually been ignored in studies of galaxy formation. The 
cosmic ray particles may originate in acceleration processes 
in high Mach number shocks in supernova remnants or could 
be produced in structure formation shock waves. We here 
focus on cosmic ray injection associated with supernovae, 
and consider only Coulomb and hadro nic interactions as lo ss 
processes for the cosmic ray particles (|Enfilin et al.ll2007l ). 

A numerical treatment of the cosmic ray component 
is rather complicated as in principle their full, in general 
anisotropic distribution function has to be modeled. Also, 
the motion of the cosmic ray fluid is tightly coupled to the 
magnetic field, which in turn is non-trivial to calculate ac- 
curately. We therefore employ the subresolution model de- 
scribed by iJubelgas et a l. (2008), which has alread y been 
successfully employed in p revious work fPfrommer et al.l 



successiully employea m p r 
l2007l . l200a : |Pfrommei|[200i ) . 

The basic assumption of our CR model l|Enfilin et al.l 
l2007l : |Jubelgas" et al.|[200^ ) is that the momentum spectrum 
of cosmic rays can be well represented by a simple power 
law of the form 



d^iV 



= Cp-"e(p-g), 



(1) 



dpdV 

where C gives the normalization, g is a low momentum cut- 
off, and Q is the power law slope. We assume that the cosmic 
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Label 


Mass resolution (gas) 


Gravitational softening 


Physics 




Ref 


5.14 X 10* M0 


0.34 kpc 


star formation, supernova feedback'" 




BH 


5.14 X 10* Mq 


0.34 kpc 


star formation, supernova feedback'" 


, AGN feedback'*) 


Wind 


5.14 X 10* Mq 


0.34 kpc 


star formation, supernova feedback'" 


, galactic winds'"' 


CR 


5.14 X 10* Mq 


0.34 kpc 


star formation, supernova feedback'" 


, AGN feedback'*"' , cosmic rays'"^' 


LowRes 


4.11 X 10^ Mq 


0.68 kpc 


star formation, supernova feedback'" 


, AGN feedback'*' 



Table 1. Overview of the simulations used in this work. T he label of each simulation w ill b e used throughout th e rest of this paper. The 
differe nt physics models we use are described in detail in (a) ISprineel fc Hernauis"3 ([200^, (b) ISprineel et al.l l l2005^ . and (c) ljubelgas et al] 
ll2008h . 



rays are dominated by protons, and that the local magnetic 
field is sufficiently tangled to effectively lock in the CRs to 
the gas. The pressure of the cosmic ray population is then 
given by 

while the number density is simply ncR ~ Cq^~°' /{a — 1). 
Here 

rn 

B„{a,b)= x''-\l-x)''-'^dx (3) 
Jo 

denotes incomplete B eta functions. Based on the findings of 
lJubelgas et al.l (|2008l ). we expect this additional pressure to 
be especially influential in low mass satellites, because it here 
can substantially reduce the density of the ISM. This should 
effectively reduce the number of low luminosity satellites. 

We model the decay of the cosmic ray population by 
accounting for coulo mb cooling and catast rophic hadronic 
losses as described in lJubelgas et all (|2008t ). Note that the 
'cosmic ray cooling' mediated by these effects can occur on 
a different timescale as the ordinary thermal cooling. In par- 
ticular, gas can end up being cosmic ray pressure supported 
after having lost much of its thermal support through ra- 
diative cooling. Finally, we adopt a simple source function 
for the injection of new cosmic ray particles, which we link 
directly to the star formation rate. This is motivated by ob - 
servations of supernova remnants (|Aharonian et aLll2006l ). 
where a large fraction of the supernovae energy is seen to 
initially appear as cosmic rays. 

2.5 Simulation set and analysis 

In Table[Tl we summarize the primary properties of the sim- 
ulations we analyze in this work. We consider four different 
high-resolution simulations of the same initial conditions, 
corresponding to the Aq-C-4 halo, but carried out with dif- 
ferent physics in the baryonic sector. Our reference calcula- 
tion (labeled 'Ref') includes star formation and supernova 
feedback as described by our multi-phase model, as well as 
ordinary radiative cooling and heating by a UV background 
that reionizes the universe by redshift z — 6. We have re- 
peated this calculation by adding in turn each of the three 
additional feedback models described above. This yields the 
three simulations 'BH' (with black hole growth and feed- 
back), 'Wind' (with the phenomenological wind model), and 
'CR' (with cosmic ray physics). Our primary simulation set 
is composed of these four simulations. They are of equal 
numerical resolution and hence allow a relatively clean as- 
sessment of the impact of the different physics components 



on the satellite population. We note that our primary aim 
in this work is not to construct a best-fitting model for the 
Milky Way, as this may require a combination of the differ- 
ent physics models and a fine-tuning of their free parame- 
ters. Rather we want to highlight the importance of different 
physics for the satellite population. 

We also briefly consider a further simulation, labeled 
'LowRes'. This is a rerun of our reference model at lower 
resolution, corresponding to 'Aq-A-5'. We use this simula- 
tion for an assessment of the numerical resolution and con- 
vergence limits of our simulations. 

To analyze the time evolution of the simulated galax- 
ies, several snapshots were stored at different times. As a 
basic analysis step, the snapshots were flrst processed by a 
group finding algorithm in order to identify individual halos. 
The group finding was done with a simple friend-of-friends 
(FOF) algorithm applied only to the dark matter particles 
with a linking length equal to 20% of their mean particle 
spacing. The gas and star particles were linked to their near- 
est dark matter particle. Next, each halo found in the first 
step was subjected to a substructure det ection procedure, for 
which we used the SUBFIND algorithm jSpringel et al.ll200ll ) 
in a version exten ded to allow a treatment of gas as well 
(|Dolag et al.ll2009l ). SUBFIND calculates the local density 
everywhere and searches for substructure candidates that 
are locally overdense. It then computes the gravitational 
potential at the positions of all particles in the candidate 
structures, and determines the subset of particles that are 
gravitationally bound. In this way, only real physical struc- 
tures are found. To avoid noise from substructures composed 
of only a few particles, only substructures containing at least 
20 particles were kept for further analysis. 

The gravitationally bound structures found by SUB- 
FIND in this way form our catalogue of simulated galajcies, 
including both central galaxies as well as genuine satellites. 
For the simulate d galaxies, we applied a stellar population 
synthesis model l|Bruzual fc Charlotll200 j ) to estimate their 
luminosities and colours, based on the formation times and 
masses of the star particles created in the simulations. We 
made no attempt to account for the metallicity dependence 
of the stellar population synthesis model or dust corrections. 



3 OBSERVATIONAL KNOWLEDGE 

Before we present our simulation results, we summarize in 
this section the most recent observational data with respect 
to the Milky Way's satellite population. We will later use 
this comprehensive catalogue of the known satellites to- 
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Table 2. Com pilati on of all presently kn own Milky Way s atell ite galaxie s. The values are taken from the following stud i es: (a ) 
Ivan den Be7gh! |l994l), (b ) Mateo ( 19981 . (d llrwin et al." ('2007h.(' d') iLiu et al.l (20081. fe) Ma rtin et al.] ll2008l^. (fllSimon fc G eha (200^ 
(s)|B ekki (.20081'). (h) iGrebellpOOO) . (il .van den Bergh (200(3,), (i'l iBelokuroTet al.. (2008). (k) iBelokurov et al.l ||2009| ). Ql lWatkins et all 
pOOgh . (ml iBelokurov et alj l201oh . All errors are from the corresponding papers. The different columns list the position in galactic 
coordinates, the proper distance, the total V-band magnitude, the surface brightness and the total estimated mass of the individual 
satellites. 



gether with predictions for their total number over the whole 
sky when comparing with our simulations. 



Table [2] gives a compilation of the properties of all 
Milky Way satellites known today. The basic data of the 
'classical' sate llites, which were already known in 1998, 
are given b v iMated (ll99S ) and were only slightly ex- 
tended using Ivan den Berghl (|2000l ) who updated the data 
on the Small and the Large Magellanic Cloud (SMC 
and LMC, respectively). Up to this time, the number 
of known satellites was just 16, but during the differ- 



ent data releasefQ (lYork et al.l I2OO0I : IWillman et~aLl I2OO2I ) 

of the Sloan Digital Sky Survey (SDSS), the number of 
known satellites increased significantly thanks to new dis- 
coveries made with the survey. Table [2] includes the new 
satellite galaxies fou n d with SDSS (W illman et al. 200 5a,i3; 
'Zuc ker et aLl l2006bl: [Belokurov et all 2006; Zucke r et al 
2006al : iBelokuroy et a,l.l|2007l: llrwin et al..,2(307.: .Walsh et al 



2001 iLiu et all l2008l: IWalsh et al.l l2008l: IBelokurov et al 
i200a . l2009l : IWatkins et al.ll2009l : IBelokurov et al.ll2010l ). us- 



|http : //www ■ sdss . org/dr6/ Index . html] 
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ing the recen t ly pu blished structural parameters given in 
iMartin et al] (|2008l '). To estimate the half light radius of 
the Large M agellanic Cloud, the formula rn = 1.68 re 
l|Martin et al.i,2008. ) was adopted. 

Up to now, there are 35 known Milky Way satellites. 
However, this sample is not complete, as the recently found 
satellites are all limited to the area of the sky covered by 
SDSS, which corresponds to a fraction of 0.194 of the full 
sky. Effectively, only this region of the whole sky ha s been 
scanned for faint satellites (see iTollerud et ahlbOOsI ). Tak- 
ing into account the det ection limits and t he sky cover- 
age of the SDSS survey, , Simon fc Gehal (|2007l ) estimate the 
number of satellite galaxies with a surfa ce brightness above 
« 28 mag/ arcsec^ (|Martin et al.l |2008| ) expected over the 
whole sky to be 57. However, more recent works favor a limit 
of 30 mag/arcsec^ jSullock et al.l 12009*). which we adopt 
throughout the rest of this paper. Using this threshold, we 
denote simulated satellite galaxies as 'observable' if their 
surface brightness exceeds 30 mag/arcsec^. 



4 ABUNDANCE OF LUMINOUS SATELLITES 

Arguably the most fundamental property of the subhalo 
population is the abundance of satellites as a function of 
luminosity. In Figure [T] we show the differential luminosity 
function of the simulated satellite galaxies for our four pri- 
mary simulations, and compare to the observations for the 
Milky Way. The latter are expressed in terms of the fitting 
formula given bv lKoposov et all (|2008l ): 

^^l0^l00.i(M^+s)^ (4) 
dMv 

The upper left panel of Figure [T] shows our result for 
the reference simulation. There is a sizable offset between 
the simulated and the observed luminosity functions, shown 
in black and green, respectively. Satellites with large stellar 
masses are even more strongly overproduced than low lu- 
minosity ones. This shows that photoheating and supernova 
feedback as included in the reference model are insufficient 
to match the observed satellite abundance. 

The other panels of the figure show the results we ob- 
tained for our alternative physics simulations. As one might 
expect, including supermassive black holes has no substan- 
tial influence on the population of satellites, yielding essen- 
tially the same result as for our reference simulation. This 
shows that environmental heating effects from quasar feed- 
back, in particular the possible quenching of star formation 
in nearby small halos, play no important role in the his- 
tory of the Aq-C halo. Interestingly, galactic outflows with 
our standard wind prescription are also unable to signifi- 
cantly improve the agreement with the observations. While 
the most luminous satellites are moderately suppressed in 
stellar mass, the satellites tend to pile up on the faint side 
of the luminosity function, yielding to a slight steepening ef- 
fect of the luminosity function, quite different from what is 
needed to match the data. Finally, the simulation including 
cosmic rays yields a substantial modification of the results. 
Here we actually obtain very good agreement with the lu- 
minosity function inferred from the observations, because 
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Figure 2. Cumulative satellite luminosity function for our dif- 
ferent simulation models, compared to the observed luminosity 
function. Note that the observed luminosity function in this plot 
includes only known satellites without any incompleteness cor- 
rection and should thus be seen as a lower limit, especially below 
—7 Mv where satellites can only be detected within the 19.4% 
sky coverage of SDSS. 



compared with the reference model the luminosity of the 
satellites is efficiently suppressed by the CR feedback. 

Another view on the satellite abundance is given in Fig- 
ure[21 where we show the cumulative abundance of the satel- 
lite population as a function of luminosity, comparing all 
four simulation results to the observations. Here, the obser- 
vational data is based only on direct observations, meaning 
that the black line should be taken as a lower limit for mag- 
nitudes lower than f» —7 My, because of the incomplete 
sky coverage of the SDSS. The figure confirms the conclu- 
sions we reached from the different results of Fig. [1] The 
'Wind' simulation is efficient in reducing the luminosity of 
large satellites of the size of the LMC/SMC, but does not 
manage to suppress the abundance of low luminosity satel- 
lites. In contrast, while the CR simulation does not reduce 
the luminosity of the brightest satellites much, it is very ef- 
ficient in suppressing star formation in low-mass subhalos, 
ultimately producing a much reduced amplitude of the lumi- 
nosity at the faint end. As a result, the CR simulation stays 
quite close to the observational data up to the completeness 
limit of the SDSS. 

Figure[2]also includes the result of the 'LowRes' simu- 
lation, shown as a green dot-dashed line. We can see that the 
population of satellite galaxies is independent of resolution 
to good accuracy up to a magnitude of « —12 Mv- Taking 
into account that the resolution limit of our high resolution 
runs is shifted by « 4 magnitudes, we expect that our sim- 
ulated satellite abundance should be numerically converged 
for satellites brighter than « —8 Mv- 

Figure [3] shows yet another way to compare the counts 
of simulated satellites with the observations. Here we use 
the maximum circular velocity of satellites, Umax, on the ab- 
scissa, because circular velocities are a good proxy for the 
(original) mass of the subhalos, but can be much more reli- 
ably measured than the mass itself. We note that such veloc- 
ity functions have already been used in the first discussions 
of the missing satellite problem, and are still frequently used 
to compare the number of observed satellites with the sub- 
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Figure 1. The differential luminosity function of the simulated satel lites (black cur v e) com pared to the observational data. The latter 
are represented by the red line, which is the fitting function given bv lKoposov et aL The green line shows the scaled luminosity 

function they obtained for the SDSS satellites while the open diamonds give an extended version using also the luminous satellites of the 
Milky Way and Andromeda. The different panels show our results for the Ref, BH, Wind and CR simulations, respectively. 



structure abundan ce in coUisionless N-body simulations (e.g. 
iMadau et al.l2008l 'l. The filled green circles in Fig.[3]show the 
raw observational data, while the open circles are a scaled 
version that accounts for the SDSS sky coverage and incom- 
pleteness. The real cumulative velocity function might thus 
be expected to lie close to the filled symbols at high circular 
velocities and to approach the open circles at low circular 
velocities. The dashed black line shows the cumulative ve- 
locity function of all satellite galaxies produced in a dark 
matter only simulation, using the same initial conditions as 
for the high resolution hydrodynamical simulations. Finally, 
the blue line in each panel shows the cumulative mass func- 
tion of all satellites containing at least 1 x 10* /i~^M0 of 
stellar mass (one star particle) in the corresponding hydro- 
dynamic simulation. 

The differences between the observations and the sim- 
ulation results for the different physics models appear large 
at first sight. This however confirms and is consistent with 
our earlier findings. In particular, the reference simulation 
and the simulation with black hole feedback overpredict the 
satellite counts for all velocities, while the wind simulation at 
least manages to give a reasonable abundance of the bright- 
est satellite systems. Again, we find the cosmic ray simula- 
tion to produce the best match to the data. Whereas there 
may still be a moderate overproduction of bright systems, 
the extrapolated faint end abundance is matched quite well. 



and, in particular, the shape of the predicted luminosity 
function is in quite good agreement with the observations. 

There is another interesting aspect of Figure|3]that con- 
cerns the comparison with the dark matter only results. 
It is a generally assumed that satellite galaxies are dark 
matter dominated. However, comparing the black dashed 
line, which shows the mass function of satellites in the cor- 
responding dark matter only simulation starting from the 
same initial conditions, with the result of the individual hy- 
drodynamic simulations, we note some sizable differences. 
The high mass satellites show clear evidence that gas cool- 
ing has led to a higher concentration of their mass profiles, 
thereby increasing their circular velocities. Despite the rel- 
atively low stellar mass content in these bright satellites, 
they hence show some structural changes due to baryonic 
effects. We note however that invoking yet stronger super- 
novae feedback may reduce these effects if the cooling rate 
is more effectively reduced. 

4.1 Kinematic results 

We close this section with an anlysis of some of the structural 
properties of the simulated satellite population. As noted 
earlier, the total mass of a satellite galaxy is difficult to mea- 
sure, so other tracers are usually used as a proxy for mass. 
An observationally readily accessible measure of this type 
is the central velocity dispersion, which is very commonly 
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Figure 3. The cumulative number of luminous satellites (with stellar mass larger than 10^ h ^Mq) as function of th eir maximum c ircula r 
velocity. We compare results for our different physics sim ulations with data for the Milky Way satellites, taken from lMadau et al] | |2008|) : 
ISimon fc Geh3 ll2007D : lMatej lll998D ^Martin et alj ||2007| ) and plotted as green filled circles together with the assumption that the circular 
velocity equals \/3 times the central velocity dispersion {Prhnack 2009). The error bars indicate a plausible range of circular velocities 
between a and 4 X cr. Open circles show the same data as the solid circles, but scaled by a factor of 5 to roughly account for SDSS 
sky coverage and incompleteness. The blue line shows the cumulative mass function of the observable satellites in each hydrodynamical 
simulation, while for comparison the black dashed line gives the mass function of all substructures in the corresponding dark matter 
only simulation. The variations in the maximum circular velocity obtained for the host halos in the different simulations are due to the 
different sizes of the stellar bulges grown in the different runs. 



used fe.g. lSimon fc Gehall2007l ). In FigurelH we compare the 
relation between central velocity dispersion and luminosity 
for our simulation satellites with data from I Simon fc Gehal 
l|2007t ). u pdated with the late st values for the known satel- 
lites from lWalker et al.l (|2009l ). 



There seems to be quite good agreement between the 
observations and the dark matter velocity dispersions of the 
Ref simulation. As the sample of measured satellites is quite 
small and has rather large error bars, the weak trend of ris- 
ing velocity dispersion with rising luminosity is not very well 
determined, but the simulation apparently follows the same 
trend. We note that the alignment of the simulated satellites 
on the left hand side of the plot is due to resolution issues 
from discreteness effects. In fact, the different 'stripes' are 
separated by just one star particle. The stellar mass in each 
stripe is therefore equal, even though some scatter in lumi- 
nosity is still present because the lu minosity was calculated 
with the iBruzual fc Chariot! (|2003l ') model, taking into ac- 
count the age and metallicity of the star particles, effectively 
giving each stellar particle its own mass-to-light ratio. 



5 HISTORY OF SATELLITE GALAXIES 

In this section we track the evolution of individual satellites, 
with the aim to study their formation paths for a range of 
individual accretion, mass loss and star formation histories. 
To this end we select nine representative satellite galaxies, 
split into groups of three that are taken from three different 
mass ranges. In Figure (5] we show our 'high mass sample', 
consisting of three satellite galaxies with a final stellar mass 
higher than 5 x 10* M©. Satellites with intermediate final 
stellar mass between 5 x 10* Mq and 5 x 10® Mq are shown 
in Figure [6l while Figure [7] gives three low-mass examples of 
satellites with a final stellar mass less than 5 x 10® Mq. 

The different panels in the three Figures [5] to [7] are or- 
ganized in the same way, and show in each case the history 
of one individual satellite (with a final stellar mass as la- 
beled in the figure). For each satellite, the top panel gives 
the redshift evolution of the dark matter, gas and stellar 
components as black, green and red solid lines, respectively. 
The middle panel shows both the star formation rate (solid 
line) and the maximum circular velocity (dashed line), as a 
function of time. Finally, the bottom panel gives the evolu- 
tion of the radial distance of the satellite to the host galaxy 
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Figure 4. Total V-band luminosity agai nst velocity dispersion . 
We show the observational data points of ISimon fc Gehai ll2007l') 
combi ned with velocity dispersions taken from IWalkereT'al] 
as black crosses, and compare to the simulated satellites 
plotted as red symbols. The two samples are in good agreement, 
except perhaps for slightly different slopes. 



(solid line) and compares this to the virial radius of the host 
{R200, dashed line). Finally, the dotted vertical line running 
through all panels highlights the epoch z — 6, which is the 
time when the UV background reionizes the universe in our 
simulations. 

The different satellite histories we have selected in 
Figs. 1 5171 show a variety of interesting evolutionary effects 
that we now discuss in turn. For definiteness, we have here 
selected the BH run, but our other simulations show qual- 
itatively very similar results. The left panel of Fig. [5] gives 
a nice illustration of the tidal and ram pressure stripping 
effects that play an important role in shaping the proper- 
ties of the satellites. It is clearly seen that the dark matter 
mass starts decreasing in distinctive steps as soon as the 
satellite has fallen into the host halo and orbits with rather 
high eccentricity. These mass stripping events correspond 
to individual pericentric passages, as is clearly seen in the 
panel that gives the distance to the host halo. Note however 
that the stellar component is not noticeably effected by this 
tidal stripping process, as expected from the fact that the 
stars of the satellite are much more concentrated than the 
dark matter. In contrast, the gas component behaves rather 
differently. Here we see clear evidence for ram pressure strip- 
ping as the dominant source of mass loss even in high mass 
satellites. Interestingly, this effect starts to set in even before 
the satellite crosses the virial radius of the host, probably 
because the gaseous halo of the host is more extended than 

-R200. 

Because ram pressure stripping depends on the density 
of the surrounding gas, one should expect to see variations 
of the mass loss rate with the radial position of the infalling 
satellite. This is indeed seen if one compares the results for 
the three different satellites shown in Fig. [5]with each other. 
The satellite on the left follows a very eccentric orbit and 
spends most of the time in the outer parts of the halo, result- 
ing in a comparatively slow gaseous mass loss. The satellite 
shown in the middle panel has an orbit with a lower ec- 
centricity that keeps it at apocentre well inside the virial 
radius, yielding a consistently higher mass loss rate. Finally, 
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Figure 8. Current stellar mass of satellites versus their maximum 
circular velocity at 2 = 6. We do not find a clear threshold in iimax 
at 2 = 6 that could be used to decide whether or not a satellite 
galaxy is able to form stars later on. Instead, we find that the 
final stellar mass shows large scatter over a considerable range of 
circular velocities at the epoch of reionization. 



the satellite shown on the right panel has a nearly circu- 
lar orbit at small radius, and loses its gas component even 
faster. 

Perhaps one of the most interesting effects seen in Fig- 
ures [5] to [7| is the effect of reionization on the star formation 
of the satellite galaxies. Quite often the simplifying assump- 
tion is made that reionization would be able to stop star 
formation in satellite galaxies entirely, yet this is clearly in 
contradiction with the findings of our simulations. In fact, all 
satellites shown in these Figures (and the same is true for the 
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Figure 5. Detailed time evolution of three representative examples of high-mass satellites, from the time of their first appearance as 
individual galaxies until 2; = 0. In each case, the top panel shows the evolution of the dark matter, gas and stellar mass components as 
black, green and red lines, respectively. The middle panel shows the star formation rate as solid line, and the maximum circular velocity 
as a dashed line (with the corresponding scale on the right j/-axis) . The bottom panel gives the radial distance of the satellite to the host 
galaxy at each timestep (solid line), and compares this to the virial radius (iJ20o) shown with a dashed line. The dotted vertical lines in 
all panels mark the 2 = 6 epoch of reionization in our simulations. 



majority of other satellites too) are producing most of their 
stars at times later than z = 6. Star formation continues 
in all examples until the gas component is removed by ram 
pressure stripping, but this time can be considerably later 
than the epoch of reionization. We hence find that the de- 
tailed orbit of a satellite galaxies tends to be more important 
for determining its final luminosity than the circular velocity 
it had at the epoch of reionization. An illustration of this can 
be seen in the histories of the satellites shown in the middle 
and right panels of Figure [G] These two satellites have quite 
similar dark matter, gas and stellar masses shortly before 
they enter the hot gaseous halo of the host galaxy, but they 
are moving on very different orbits. The satellite in the mid- 
dle panel is on a relatively circular orbit, resulting in a small 
effect of tidal stripping on the dark matter component and 
no noticeable effect on the stellar component. In contrast, 
the satellite shown in the right panel is on a highly eccentric 
orbit with e > 10. This satellite dives deeply into the host 
halo, resulting in a tidal radius that is even smaller than 
the characteristic radius of the stellar component. Because 
of this, the stellar component looses nearly 90 % of its mass 
due to tidal effects. 



In Figure [S] we compare the final stellar mass of the 
satellites against the maximum circular velocity «max they 
had at redshift z — 6. This provides another way to test the 
popular hypothesis that the mass at the time of reioniza- 
tion determines the final luminosity of a satellite galaxy. We 
show results for the simulations Ref, BH and CRQ. While 
most satellites with circular velocities below ~ 20 kms~^ are 
strongly suppressed in stellar mass, there are a few objects 
with such low circular velocities that have stellar masses as 
high as 10* M©, or even 10** M©, at the present epoch. In 
the range of ~ 20kms~^ to ~ 30kms~^, no sharp cut-off 
is readily apparent that could be identified with reioniza- 
tion. Instead there is considerable scatter in the relation be- 
tween final stellar mass and maximum circular velocity at 
z — 6. We note that the simulation with cosmic rays shows a 
strong suppression in the stellar mass for low circular veloci- 
ties when compared with the other simulations, as expected 
from our luminosity function results. 

t Unfortunately, all snapshots before 2 = 2.7 of the simulation 
Wind where accidentally deleted, making this comparison impos- 
sible for this model. 



© 0000 RAS, MNRAS 000, 000-000 



12 M. Wadepuhl & V. Springel 



1 0.5 




10" 
10" 
10= 
10" 
10' 
lO** 
10' 

o'^ 

0.20 
0.15 
0.10 
0.05 

%88 

500 
400 



S 300 

DC 

200 



10 



5 3 2 1 0.5 



100 







- 


' ' ' — ' 





III — ■ 





100 
80 

60 ^ 
40 . 
20 




0.1 



1.0 




Figure 6. The same as Fig. (5] but for three intermediate mass satelHtes. 



A complementary of view of the above relation is shown 
in Figure!^ where we plot the current maximum circular ve- 
locity of satellites against their current stellar mass. Differ- 
ent symbols encode the circular velocity they had at redshift 
z = 6, where satellites with Wmax > 16 kms~^ at z = 6 are 
shown as red diamonds while satellites below this threshold 
are shown as green stars. Note that there is considerable 
overlap between the regions occupied by the different sym- 
bols, showing again that the correlation between the circular 
velocity at the epoch of reionization and the current stellar 
mass is not overly strong. In particular, there are some ex- 
amples (especially in the BH simulation) of satellite galaxies 
with very low circular velocity at 2: = 6 which were never- 
theless able to form many stars later on and to turn into 
reasonably luminous satellites. 

Finally, in Figure [TOl we directly compare the results of 
the different physics simulations with each other, in terms 
of the relationship between Wmax and the stellar mass at the 
time of reionization. Interestingly, we see that AGN feed- 
back is in fact able to reduce the stellar mass formed at 
high redshift in many of the progenitor systems of today's 
satellites, even though the effect is considerably weaker than 
for cosmic rays. As we have already seen, the influence of 
BH feedback tends to become however low at later times, so 
that the present day properties of satellites are only mildly 
effected. This is presumably because most satellites do sim- 
ply not grow a massive black hole, but they are nevertheless 



affected at high redshift by the seed black hole that is in- 
jected into their halo. 

In Figure llll we plot the baryon fraction of satellites 
against Umax, at the present epoch. The baryon fraction is 
here simply defined as the total bound baryonic mass rela- 
tive to the total bound mass of a halo. It is interesting to 
compare this value with the universal cosmic baryon fraction 
Qb / {i^b + ^dui) = 0.16 (shown as dashed horizontal line). As 
one expects, the baryon fraction is usually lower than the 
cosmic baryon fraction, especially for very low mass satel- 
lites that have lost most of their gas and did not form many 
stars either. However, in the simulation with comparatively 
weak feedback, some satellites have also baryon fractions 
above the cosmic mean value. These are satellites which lost 
a lot of dark matter through tidal stripping whereas they 
could hold on to most of their stars. Both the Wind and CR 
models are leading to considerably reduced baryon fractions 
in low mass satellites. In the former case, this is readily ex- 
pected as a signature of the winds. In the latter, it is because 
more baryons stay in a diffuse gaseous phase, allowing them 
to be more easily ram-pressure stripped. 

An analysis of the evolution of the baryon fraction be- 
tween z — 6 and 2 = is given in Figure 1121 We here 
only show results for the Ref simulation, as the qualita- 
tive behavior of the other simulations is similar. We use two 
symbols for each satellite, one showing the data point at 
z — 6 (red stars), while the corresponding values at 2; = 
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Figure 7. The same as Fig. [5] but for three low mass satellites. 



are given by green stars. Every pair of points belonging to 
the same satellite is connected by a dotted line. Most satel- 
lites with circular velocities below ~ 20kms^^ at 2 = 6 
lower their baryon fraction substantially until the present 
epoch, and they also do not tend to grow much. In contrast, 
most larger satellites tend to keep their baryon fraction or 
increase it slightly, often accompanied by a significant in- 
crease in Umax- In the most massive satellites, part of this 
increase stems from modifications of the inner rotation curve 
due to the formation of a quite concentrated stellar compo- 
nent, i.e. these satellites are not really bona-fide dark matter 
dominated systems as often assumed. We note however that 
the threshold at ~ 20kms~^ is not sharp; there are still 
many examples of satellites with an initially high «max that 
end up as low mass satellites with a stripped baryonic com- 
ponent. 

The last quantity we analyze in this section are the cu- 
mulative star formation histories of our satellites, as shown 
in Figure 1131 The solid black line shows the total cumula- 
tive star formation history of all satellites in the final virial 
radius, normalized by their total final stellar mass. The gray 
shaded area gives the la scatter around this mean for the 
ensemble of all satellite star formation histories. The ver- 
tical dotted, dashed and dot dashed lines mark the times 
when 10 %, 50 % and 90 % of the stars present at z = were 
formed. Finally, the dashed blue line repeats the result of 
the Ref simulation in all the panels corresponding to the 



other simulations, in order to ease a comparison between 
them. As we have already seen in the other results, AGN 
feedback shows little effect on the cumulative star forma- 
tion history of the satellites. The Wind model on the other 
hand leads on an earlier production of the bulk of the stars, 
which is what one would expect if galactic outflows are ef- 
ficiently removing gas from star-forming dwarf galaxies and 
are thus shutting down star formation earlier. In contrast, 
the CR model shows exactly the opposite effect. Due to 
the additional pressure component, the galactic gas has a 
lower overall cooling rate. This hampers star formation in 
low mass systems but does not by itself remove significant 
amounts of fuel for star formation; the latter can however 
be achieved by ram pressure stripping. Thus, star formation 
shifts to considerably later times in the CR run than in any 
of the other models. 

Interestingly, the scatter around the mean history of 
the satellites is also modified by the different physics. The 
Wind simulation shows a rather small scatter, presumably 
because most satellites form their stars in the first significant 
phase of star formation at high redshift, which is terminated 
quickly and for the most part coevally. In the case of the CR 
simulation, much of the gas is not removed by the primary 
feedback process itself, but instead is affected by stripping 
processes at intermediate and low redshifts, after the satel- 
lites have fallen into the parent halo. This means that the 
individual infall history of each satellite is of larger impor- 
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Figure 9. Current stellar mass versus current maximum circular 
velocity. Satellite galaxies with i)max > 16 km/s at z = 6 are 
shown as red diamonds, and satellites below this limit as green 
stars. There is again no evidence that this threshold or one nearby 
is able to distinguish between satellites with or without a signifi- 
cant amount of stars. 



tance in this model, leading to a higher overall variability in 
the star formation history of the satellites. 



6 SCALING RELATIONS 

In this section, we investigate in more detail how the prop- 
erties of our simulated satellites scale with their size. Where 
possible, we compare with observational results and other 
theoretical predictions. We want to caution however that 
especially our smallest luminous satellites are pretty close 



Figure 10. Stellar mass versus maximum circular velocity at 
z = 6, compared for different simulation models. Interestingly, 
a weak influence of AGN feedback on the satellite population is 
found at high redshift, but this diflference largely vanishes later 
on, as most of the satellites are simply not able to grow a large 
supermassive black hole. In contrast, the additional pressure com- 
ponent of cosmic rays affects all satellites more strongly, and here 
the effect remains large in small galaxies even down to 2 = 0. 



to our resolution limit. While the satellites above the detec- 
tion limits of SDSS should be sufficiently well resolved in 
our high resolution simulation to give reliable results, a con- 
siderable numerical uncertainty persists, a fact that should 
be taken into account in interpreting the results. 

We begin with the scaling relations derived by 
IWoo et al.l (|2008h for local group dwarf galaxies. We focus 
on the relations between stellar mass and circular velocity, 
and stellar mass and star formation rate, as they have the 
highest statistical significance and are thus best suited to 
benchmark the simulation results. Figure fT^ shows these two 
relations i n separate panels, comparing in each case the fits 
derived by I Woo et all (|2008l ') with our simulation data. The 
correlation between stellar mass and maximum c ircula r ve- 
locity is comparatively tight, in fact, IWoo et al.l (|2008l ') cite 
a correlation coefficient of 0.94 for the observations, which 
are represented by the solid red line. The simulated satel- 
lites show a similarly strong correlation (formally yielding a 
correlation coefficient of 0.94), but the results for the Ref 
simulation are slightly offset towards higher stellar masses. 
The CR results (shown with magenta symbols) agree con- 
siderably better. Our lowest mas s satel lites start to deviate 
from the fit given bv lWooet al.l Hooi) but the differences 
are of comparable size as in some observed systems such 
as Ursua Minor, and the growing scatter in the numerical 
results also indicates that resolution effects start to play a 
role. 

The right panel in Fig. [TJ] compares the correlations 
between stellar mass and current star formation rate. Here 
we find a much worse ag reement with the observational re- 
sults of I Woo et all (|2008l ). which are again characterized by 
a remarkably good correlation (with coefficient 0.96). While 
our results bracket the observationally inferred relation, the 
scatter is large and the formal correlation is only 0.72. In 
addition, most of our satellite galaxies show only a vanish- 
ingly small star formation rate at the present epoch, and 
those systems were omitted from the plot. But again, the 
CR results seem to agree better with observations. Simu- 
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lations with better mass resolution will be needed to shine 
mor e light on this p otential discrepancy. 

ISimon fc Gehal (|2007) calculated the mass-to-light ra- 
tio for the sample of Milky Way satellites found in the 
SDSS, obtaining values ranging from about 100 to 1200, 
with a mean of ~ 380. Doing the same calculation for 
the whole sample of known satellites resulted in values be- 
tween 1.5 and 1200 with a mean of ~ 170. This suggests 
that the faint satellites discovered with the SDSS are even 
more dark matter dominated than the more luminous 'clas- 
sical' satellites. For the complete sample of simulated lu- 
minous satellites, we obtain mass-to-light ratios between 
12 (11, 16, 10) and 13000 (18000, 23000, 18000) with a mean 
of ~ 1500 (1350, 1550,3500) for the Ref (BH, Wind, CR) 
simulation. As mentioned earlier, very small satellites are 
strongly affected by numerical effects and thus the very large 
mass-to-light ratios we find for these satellites may be un- 
reliable. Also, the mean value may be biased high by the 
large number of small satellites. If we restrict the mass-to- 
light ratio calculation to satellites with Mv < —8.0, we ob- 
tain a mean of 232 (217,315,345), which is much closer to 
the observed values. If we select only satellites with a sur- 
face brightness fj, < 30, then the mean mass-to-light ratio 
is 33 (33,73,33), which is about 5 times smaller than the 
observed mean value for this selection. 

This difference in the mean mass-to-light ratio can also 
be seen from the "Mateo Plot" shown in Figure 1151 which 
compares the luminosities of the satellite galaxies with their 
mass-to-light ratio in units of the solar mass-to-ligh t ratio. 
In the original paper where this plot was introduced, iMated 
(ll99Sr) overplotted the function 

i°S(]^ = 2-5 + ioV(i^/^«)' (5) 

which we also included as the dark red dot-dashed line in 
Fig. 1151 To guide the eye, we also simply scaled this function 
by a factor of « 5.2 and plotted it again as the purple dot- 
dashed line. It can be seen that this scaled function fits the 
simulated galaxies of the Ref simulation very well, while the 
observed systems (shown with red and green triangles) are 
well described by the original function. This result is consi- 
tently reproduced by all simulations. Only the constant hor- 
izontal offset between observations and simulations changes 
by ~ 10% while the general trend remains the same. 

We note that the nearly constant offset between the 
simulated and observed satellite galaxies could in part be 
caused by the rather uncertain procedure applied to esti- 
mate the total mass of observed satellites. This effectively 
involves an extrapolation to the outer edge of the satellite, 
which is uncertain. An alternative would be that the simu- 
lated galajdes simply contain fewer stars than expected for 
an observed satellite of the same mass. However, the Ref 
simulation already has comparatively weak feedback, and 
allowing for brighter satellites by a constant factor would 
cause the most luminous satellites, which are in good agree- 
ment with the Magellanic Clouds, to become too bright. 
Furthermore, making the star formation more efficient in all 
satellites would shift the points in Fig. [TJ] both down and 
to the right, hence spoiling the good agreement with the 
location of the break in the observed relation. 

In Figure fTUl we show the relation between photomet- 
ric surface brightness of all simulated dwarf galaxies inside 



1 .0000 r 




0.0001 I I 

10 20 30 40 50 

Figure 12. Evolution of the baryon fraction and circular velocity 
between z = 6 and z = 0. For each satellite, we mark the high 
redshift z = 6 data with a red star and the z = value with a 
green star, and we connect each pair of two points with a dotted 
line. There is clearly a pivotal maximum circular velocity of about 
~ 25kms~^ below which most satellites lose a large fraction of 
their baryons, while they retain their baryon fraction above this 
threshold. 

a sphere of radius 350 kpc centred on a fiducial position of 
the Sun. The Sun was assumed to lie 8.5 kpc away from 
the centre of the galaxy, in the central plane of the stellar 
disc. As can be seen from the relatively large scatter of the 
plot, the simulation produces also satellites that are well 
above the SDSS surface brightness detection limit. Count- 
ing the galaxies with a photometric magnitude brighter then 
30 mag/arcsec^ (dashed red line) and a distance smaller 
than 280 kpc (dot dashed blue line) results in observable 46 
(77, 18, 70) satellites for the Ref (BH, CR, Wind) sim- 
ulation. This is actually in reasonable agreement with the 
prediction of 57 satellites for the Milky Way. The cutoff ra- 
dius of 280 kpc was chosen as a compromise between the 
measured distances to all known satellites, which reach up 
to « 1 Mpc, and the virial radius of r2oo = 238 kpc of the 
simulated host galaxy. This is also th e same cut off radius 
that has been use d in previous work (|Koposov et al.l [20081 : 
iMaccid et alj2O10l ). although we note that some studies have 
adopted a different choice (e.g. _Diemand ct al. 2007). The 
rather small number of satellites classified as 'observable' 
for the CR simulation can easily be explained by the lower 
luminosity function shown in figure [1] 

Finally, we consider the relation between the dark mat- 
ter masses of our simulated satellites with their stellar mass 
and luminosity, as shown in Figure ITTl For each satellite, we 
plot the dark matter mass with different symbols, both at 
the epoch of accretion and at the present epoch. T o simplify 
a comparison with figure 5 of iMaccid et al.l ()2010l ). we used 
exactly the sa me axis range in our plot as they did. Unlike 
in the results of lMaccid et al.l(|2010h . we find a clear bend in 
the relation, meaning that our satellites tend to have higher 
stellar masses, especially at the low mass end, than the satel- 
lites of lMaccio et al.l l|2010f ) . The latter results are based on a 
semi-analytic model where the orbits of an infalling satellite 
are estimated based on a random choice of plausible infall 
parameters. It is possible that this explains the discrepancy, 
or that it originates in approximate treatments of tidal or 
ram pressure stripping in the semi-analytic model. In future 
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Figure 11. Baryon fraction of 2: = satellite galaxies as a function of their circular velocity. We show results for our four primary 
simulation models, and define the baryon fraction in terms of the bound particles identified by SUBFIND for each satellite. The horizontal 
dashed line gives the universal cosmic baryon fraction. 



work, it will be interesting to inter-compare direct hydro- 
dynamical simulations and the semi-analytic models on a 
satellite by satellite basis, in order to better understands 
the origin of these differences in the predictions. 



7 CONCLUSIONS 

In this work, we studied a set of high-resolution hydrody- 
namical simulations of the formation of a Milky Way sized 
gala:xy, starting from cosmological initial conditions. Such 
simulations are now able to reach sufficiently high resolution 
to directly resolve the formation of the small dwarf galaxies 
that orbit in the halo, thereby allowing studies of the missing 
satellite problem and of the properties predicted by simula- 
tions for the population of satellite galaxies. These galax- 
ies are especially interesting both because the dark matter 
substructure abundance is a fundamental challenge for the 
ACDM cosmology, and because the low star formation ef- 
ficiencies of the satellites provide crucial information about 
the physics of feedback. 

We have therefore repeated our simulations using dif- 
ferent models for feedback physics, with the goal to test 
the sensitive of the results for the satellites with respect to 
these physics assumptions. In the Ref model, we consid- 
ered only star formation and SN feedback, together with 
instantaneous reionization at 2 = 6. The three other mod- 
els included additional processes like AGN feedback (BH), 



wind driven galactic outflows (Wind) and the generation 
and decay of cosmic rays (CR). Not unexpectedly, the BH 
model showed no significant differences compared to the ref- 
erence Ref model, as most of the satellites are simply too 
small to grow a large supermassive black hole and are rarely 
affected by strong quasar feedback in neighbouring galaxies. 
In contrast, the Wind model showed a significant reduction 
of the number of high mass satellites, but did not give a sig- 
nificantly different abundance of low mass systems. The CR 
model had exactly the opposite effect as it did not change 
the high mass satellites but suppressed star formation in 
low mass satellites. This made the cosmic ray model most 
successful in matching the faint-end of the observed satel- 
lite luminosity function. Our results further suggest that a 
combination of the Wind and CR feedback models should 
be able to yield a nearly perfect match of the luminosity 
function. 

The total number of satellites observable with an SDSS- 
like survey covering the w hole sky has been estimated to be 
57 l|Simon fc GehallioOTi ). Interestingly, imposing the same 
surface brightness detection threshold on all of our simu- 
lated systems yields a prediction of 77 observable satellites 
for our BH model, which is only moderately higher than 
the observations despite the fact that this simulation over- 
predicts the satellite luminosity function considerably. For 
our CR model instead, the number drops considerably, to 
18, perhaps caused in part by an overprediction of the ef- 
fective stellar radii of the satellites, which could easily arise 
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Figure 13. Average cumulative star formation histories of all satellite galaxies, for all four primary simulation models. In each panel, 
the solid line shows the average star formation history, with the grey bands mark the 1 rr scatter of the distribution. The vertical dashed 
and dotted lines give the times when 10%, 50%, and 90% of all the stars have formed. To ease the comparison between the different 
simulations, the result of the Ref simulation is repeated in all the panels as a dashed blue line. 
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Figure 14. The left panel shows the relation between stellar mass and maximum circular velocity at the presen t time , while the right 
panel gives the relation between stellar mass and star formation rate. The red lines give the best fits IWoo et al.l 1 I2OO8I) derived for the 
observational data. We here included all simul ated dwarf galaxi es within the full high resolution region of the Ref simulation (in black) 
and of the CR simulation (in magenta), since IWoo et al.] 1I2OO8I ) did also include dwarfs outside of the virial radius of the Milky Way. 
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Figure 15. The relation bet ween mass-to -light ratio and lumi- 
nosity of satellite galaxies (see lMatedll998h . The observed satel- 
lites are plotted as triangles, with red symbols marking satellites 
known before SDSS and green symbols marking the newly discov- 
ered satellites discovered in the SDSS data. Diamonds give our 
simulated satellites (Rep model), colour-coded as magenta if their 
surface brightness is high {/i > 30) or as dark blue if it is low. 
The dot-dash ed lines represent the fitting function suggested by 
iMateol l|l99^, for the observed sample (red) and shifted upwards 
by a factor of 5 (dark purple) to match the simulated sample. 
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Figure 16. Photometric V-band surface brightness of the simu- 
lated satellites (black crosses), inside a sphere of radius 350 kpc 
around a fiducial position of the Sun in the simulations. The red 
line shows the surface brightness detection limit of the SDSS sur- 
vey, while the blue vertical line gives the radial cut at 280 kpc 
that we frequently used in this work. There are 46 Satellites be- 
low this limit within a radius of 280 kpc which is relatively close 
to the 57 satellites predicted for an all-sky extrapolation of the 
observational data. The green stars give the observed satellites 
with well determined surface brightnesses. 
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Figure 17. Relation between dark matter satellite mass and stel- 
lar mass, or stellar luminosity, respectively. For each satellite we 
show both the dark matter mass today and at redshift z = 6. Thi s 
plot may be compared directly to figure 5 of lMaccio et al 
which is representative for semi-analytic models constructed to 
describe the satellite population. Unlike in their results, there is 
clearly some curvature found with the hydrodynamical simulation 
in this relation. 



from the limited spatial resolution of our simulations. In any 
case, this stresses that a large number of additional satellites 
may actually still be hidden just below the surface bright- 
ness limit of the SDSS (see also iBuUock et al.|[2009l ). 

Our simulations have also highlighted the relative im- 
portance of some of the evolutionary aspects of satellite 
galaxies. In particular, we do not find a very distinctive 
mark of the epoch of reionization on the satellites, and most 
satellites continue their star formation activity in our sim- 
ulations to much lower redshift than z — 6. This suggests 
that simplified treatments of sateUite histories, where rela- 
tively high cooling thresholds due to a ionizing UV back- 
ground are invoked, are not particularly realistic. Our sim- 
ul ation results agree m uch better with the scenario outlined 
m IStrigari et ai](|2007l ). which in fact resembles many of our 
simulation findings quite closely. 

We find that the observed relationship between V-band 
luminosity and velocity dispersion is quite well reproduced 
by our simulations, albeit with large scatter. The small 
amount of reliable observational data for the velocity disper- 
sions leaves it unclear at present whether the larger scatter 
we find indicates a problem of the simulations or whether is 
is also present in reality. What is comparatively clear though 
is that the observed relation between stellar mass and max- 
imum circular velocity is really tight, a finding that is also 
reproduced by our simulation results. On the other hand, the 
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correlation between present-day stellar mass and star forma- 
tion rate seen in our simulations seems to be not nearly as 
well-defined as in the observational data. This is related to 
the fact that we do not find a good correlation between the 
present stellar and gaseous masses; many simulated satellites 
have comparable stellar masses but differ in their gas frac- 
tions by huge factors. Gas-rich and completely gas-depleted 
satellites coexist in the same total and stellar mass regime, 
rendering a tight correlation with the star formation rate 
unlikely. 

But perhaps the most significant discrepancy between 
the simulation results and observations lies in the inferred 
mass-to-light ratios. The mass-to-light ratios of the simu- 
lated galaxies are off by about a factor of 5 when compared 
at face value to the observational estimates. This means that 
they are either too massive, or too faint for their mass. The 
discrepancy could also be caused by a systematic underes- 
timate of the total satellite masses in the observations. Due 
to the difficulty of reliably determining the 'outer edge' of 
the dark matter halo of an orbiting satellite, this possibility 
cannot be easily excluded. 

In summary, we find that the current generation of cos- 
mological hydrodynamic simulations is able to explain many 
properties of the observed satellite population surprisingly 
well. We have shown that different feedback physics affects 
the satellite population strongly, with respect to quantities 
such as luminosity function, scaling relations, or star for- 
mation histories. This emphasizes the significant potential 
of "near-field cosmology" within our Local Group to inform 
the general theory of galaxy formation. Our work has also 
shown that it is not necessarily the physics of cosmic rcion- 
ization and supernova feedback alone that is responsible for 
resolving the missing satellite problem. In fact, the role of 
reionization has probably been grossly overstated in many 
previous works, while other important feedback, such as cos- 
mic rays, has been ignored. It will therefore be very inter- 
esting to refine the hydrodynamical simulations further in 
future work, and to make them more faithful in capturing 
all the relevant physics. 
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